#pragma once

#include <iostream>
#include <vector>
#include "Polynomial.h"
#include <eigen3/Eigen/Dense>
#include <algorithm>

using namespace std;
using namespace Eigen;

class Bspline_base
{
    protected:
    vector<double> t;
    int k;
    public:
    Bspline_base(){};
    Bspline_base(vector<double> _t, int _k)
    {
        if (_k < 0)
        {
            cout << "error! k need >= 0. " << endl;
        }
        if (_t.size() - 2 != _k)
        {
            cout << "error! length of t not fit." << endl;
        }
        t = _t;
        k = _k;
    }
    Bspline_base(double _t_i, int _k)    //cardinal 
    {
        if (_k < 0)
        {
            cout << "error! k need >= 0. " << endl;
        }
        k = _k;
        for (int i = 0; i <= k+1; i++)
        {
            t.push_back(_t_i - 1 + i);
        }
    }
    double operator()(double x)
    {
        if (k == 0)
        {
            if (x > t[0] && x <= t[1]) return 1;
            else return 0;
        }
        vector<double> t1(t.begin(),t.end()-1);
        vector<double> t2(t.begin()+1,t.end());
        Bspline_base B1(t1,k-1);
        Bspline_base B2(t2,k-1);
        return B1(x) * (x - t[0]) / (t[k] - t[0]) + B2(x) * (t[k+1] - x) / (t[k+1] - t[1]);
    }
};

class Bspline
{
    private:
    vector<double> x;
    vector<double> f;
    int N;
    vector<double> coef;    //a[i]
    vector<Bspline_base> B; //B  

    public:
    double operator()(double x)
    {
        if (coef.size() == 0) cout << "请先运行cubic_inter()" << endl;
        double ret = 0;
        for (int i = 0; i < N+2; i++)
        {
            ret += coef[i] * B[i](x);
        }
        return ret;
    }

    public:
    Bspline(){};
    Bspline(vector<double> _x, vector<double> _f)
    {
        N = _x.size();
        if (N != _f.size()) cout << "Error! 点与值不等长。" << endl;
        x = _x;
        f = _f;
    }
    void cubic_inter(double dfa,double dfb);
    void quadratic_inter(double fa,double fb);

};

void Bspline::cubic_inter(double dfa,double dfb)
{
    Eigen::MatrixXd M = Eigen::MatrixXd::Zero(N, N);
    Eigen::VectorXd rhs = Eigen::VectorXd::Zero(N);

    M(0,0) = 2;
    M(0,1) = 1;
    rhs[0] = 3 * f[0] + dfa;
    M(N-1,N-2) = 1;
    M(N-1,N-1) = 2;
    rhs[N-1] = 3 * f[N-1] - dfb;

    for (int i = 1; i < N - 1; i++)
    {
        M(i, i) = 4;
        M(i, i - 1) = 1;
        M(i, i + 1) = 1;
        rhs[i] = 6 * f[i] ;
    }

    Eigen::VectorXd a = M.lu().solve(rhs); 

    coef.resize(N+2);
    coef[0] = a[1] - 2 * dfa;
    for (int i = 0; i < N; i++)
    {
        coef[i+1] = a[i];
    }
    coef[N+1] = a[N-2] + 2 * dfb;

    for (int i = -1; i <= N; i++)
    {
        Bspline_base B_i(x[0]-1+i,3);   //发生偏移的BUG
        B.push_back(B_i);
    }

}

void Bspline::quadratic_inter(double fa,double fb)
{
    Eigen::MatrixXd M = Eigen::MatrixXd::Zero(N, N);
    Eigen::VectorXd rhs = Eigen::VectorXd::Zero(N);

    M(0,0) = 5;
    M(0,1) = 1;
    rhs[0] = 8 * f[0] - 2 * fa;
    M(N-1,N-2) = 1;
    M(N-1,N-1) = 5;
    rhs[N-1] = 8 * f[N-1] - 2 * fb;

    for (int i = 1; i < N - 1; i++)
    {
        M(i, i) = 6;
        M(i, i - 1) = 1;
        M(i, i + 1) = 1;
        rhs[i] = 8 * f[i] ;
    }

    Eigen::VectorXd a = M.lu().solve(rhs); 

    coef.resize(N+2);
    coef[0] = 2 * fa - a[0];
    for (int i = 0; i < N; i++)
    {
        coef[i+1] = a[i];
    }
    coef[N+1] = 2 * fb - a[N-1];

    for (int i = -1; i <= N; i++)
    {
        Bspline_base B_i(x[0]-0.5+i,2);   
        B.push_back(B_i);
    }
}
